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ABSTRACT 

The dynamical friction experienced by a body moving in a gaseous medium is different 
from the friction in the case of a collisionless stellar system. Here we consider the orbital 
evolution of a gravitational perturber inside a gaseous sphere using three-dimensional 
simulations, ignoring however self-gravity. The results are analysed in terms of a 'local' 
formula with the associated Coulomb logarithm taken as a free parameter. For forced 
circular orbits, the asymptotic value of the component of the drag force in the direction 
of the velocity is a slowly varying function of the Mach number in the range 1.0-1.6. 
The dynamical friction timescale for free decay orbits is typically only half as long as 
in the case of a collisionless background, which is in agreement with E.C. Ostriker's 
recent analytic result. The orbital decay rate is rather insensitive to the past history 
of the perturber. It is shown that, similar to the case of stellar systems, orbits are not 
subject to any significant circularization. However, the dynamical friction timescales 
are found to increase with increasing orbital eccentricity for the Plummer model, whilst 
no strong dependence on the initial eccentricity is found for the isothermal sphere. 

Key words: galaxies: clusters: general - galaxies: kinematics and dynamics - galaxies: 
star clusters - hydrodynamics - waves. 



1 INTRODUCTION 



A gravitating body moving through a background medium 
suffers from a drag force due to the interaction with its 
own induced wake. Generally speaking, the medium may 
be composed of collisionless matter, and/or gas. For colli- 
sionless backgrounds, the classical Chandrasekhar formula 
has proved useful in determining the orbital decay of galac- 
tic satellites and globular clusters around spherical systems 
such as elliptical galaxies (e.g. Lin & Tremaine 1983; Cora, 
Muzzio & Vergne 1997, and references therein), even though 
it was inferred for uniform media. For a uniform and infinite 
gaseous medium, the gravitational drag on a body moving 
at constant velocity on a straight-line orbit, has been es- 
timated both by linear theory (e.g. Ruderman & Spiegel 
1971; Rephaeli & Salpeter 1980; Ostriker 1999) and from 
numerical experiments in different settings (e.g. Shima et 
al. 1985; Shima et al. 1986; Shankar, Kley & Burkert 1993; 
Kley, Shankar & Burkert 1995; Ruffert 1996). All these au- 
thors consider a uniform medium because they are mainly 
interested in the accretion flow past the body. Here we will 
investigate the gaseous drag and the sinking decay of a light 



body orbiting on a gaseous sphere which is initially in hy- 
drostatic equilibrium with a given gravitational potential. 
In fact, dynamical friction for gaseous and spherical back- 
grounds has its relevance to many astrophysical studies, e.g. 
galaxies in galaxy clusters, the relaxation of young stellar 
clusters, or the dynamics of massive black holes within their 
host galaxies. It will be useful to check in which conditions 
the drag formula obtained in linear theory for uniform me- 
dia works for spherical backgrounds, in analogy to stellar 
systems. 

The question of dynamical friction for gaseous and 
spherical backgrounds has received new interest in connec- 
tion with recent studies suggesting that Galactic dark mat- 
ter may be in the form of cold molecular clouds either dis- 
tributed in a disc (Pfenniger, Combes & Martinet 1994) or 
in a quasi-spherical halo (e.g. De Paolis et al. 1995; Ger- 
hard & Silk 1996; Walker & Wardle 1998; Walker 1999; De 
Paolis et al. 1999; Sciama 2000). For certain parameters of 
these clouds, the assumption of collisionless matter is no 
longer valid. Other authors have proposed, based on obser- 
vational astrophysical grounds, that dark matter may be 



self-interacting with a large scattering cross-section ( Spergel 



Steinhardt 2000). It is therefore important to find out 
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the dynamical implications of having spherical halos of col- 
lisional matter and, in particular, its effect on the decay of 
satellite galaxy orbits. 

Deeper physical insight into the question of how dy- 
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namical friction is affected by collisions may be provided 
by con sidering; continuous gaseous media. Recently, Ostriker 
(1999) pointed out that for supersonically moving bodies 
the drag in a uniform gaseous medium is more efficient than 
in the case of collisionless media that are described by the 
standard Chandrasekhar formula; and is non-vanishing for 
subsonic bodies. The same feature is observed in the case 



of a perfectly absorbing body (Ruffert 1996). As a conse- 



quence, satellite galaxies may experience more rapid decay 
in halos made up of molecular clouds, especially in earlier 
epochs before most of the gas has turned into stars. 

In the standard model of cosmological structure forma- 
tion, forming substructures continuously exchange orbital 
energy with both collisionless material and with their sur- 
rounding hot medium due to dynamical friction. An under- 
standing of the processes that can produce velocity and spa- 
tial bias where the concentration and velocity dispersion is 
different for different galactic populations, e.g. mass segrega- 
tion in clusters of galaxies, is essential for the interpretation 
of observational data as well as for semi-analytical models 
of the formation of galaxies and galaxy clusters. Clumps of 
cold gas moving through a hotter medium may suffer a hy- 
drodynamical drag that is stronger than the gravitational 
drag (e.g. Tittley et al. 1999). In this work, however, we 
focus on a purely gravitational perturber. 

The linear response of a spherical system to an orbiting 
gravit ational perturber was carried out by Balbus & Soker 
(1990) motivated by the problem of excitation of internal 
gravity waves by galaxies in the core of clusters. In this 
paper we concentrate on the dynamical friction of a rigid 
perturber orbiting in a gaseous medium that shows a spher- 
ical density distribution with a concentration towards the 
centre. Initially the gas is at rest and in hydrostatic equilib- 
rium with an external spherical gravitational potential. The 
gravitational perturber is assumed to be very large (typically 
the softening gravitational radius is taken a few times the 
accretion radius) and unable to accrete matter. We do not 
model the existence of a physical surface on the body, i.e. no 
inner boundary condition on the body has been imposed. A 
discussion on the importance of accretion is given in Section 
2. Our aim is to compute the evolution and sinking time of 
the perturber towards the center. 

The assumptions of infinity and homogeneity of the 
medium, that was required in Ostriker's analysis, are re- 
laxed, and, in addition, the contribution of non-local effects 
on the drag force are taken into account. Also, an additional 
difference is that we deal with extended perturbers whose 
radius is a significant fraction of the size of the background 
system. Since the pertu rbers in most cases of inter est do not 
present circular orbits ( van den Bosch et al. 199£ ) , and also 
motivated by the fact that cosmological simulations show 
that the majority of satellite orbits have quite large eccen- 
tricities, we consider the dependence of the sinking times on 
the eccentricity of the orbits. 

There are some similarities with the dynamical evo- 
lution of binary star cores embedded in a common enve- 
lope (e.g. Taam, Bodenheimer & Rozyczka 1994; and refer- 
ences therein). In those systems, however, the gravitational 
torques of the two cores are able to spin up the envelope 
to a high rotation (see also Sandquist et al. 1998). The fact 
that the masses of the stars and the envelope are compa- 
rable is crucial for the subsequent evolution of the system. 



Here we explore a rather different case in which there is a 
single rotating body and its mass is significantly less than 
that in the gaseous background mass. Recall that neither 
wind is injected in the gaseous component nor radiation is 
considered in the present work. 

The paper is organized as follows. In Section 2, we high- 
light the differences between dynamical friction in collision- 
less and collisional media; the importance of having mass 
accretion by the perturber is also addressed. The descrip- 
tion of the model is given in Section 3. In Section 4 the 
orbital decay rate is computed numerically for a body orbit- 
ing around a gaseous sphere, and is compared with a 'local' 
estimate of the drag force. Implications and conclusions are 
given in Section 5. 



2 THE EFFECTS OF COLLISIONS AND 

ACCRETION ON DYNAMICAL FRICTION 

2.1 Comparing the drag force in collisionless and 
collisional backgrounds: a statistical discussion 

Dynamical friction may be pictured as the response of the 
system which tends to equipartition. The most classical 
problem is the motion of a 'macroscopic' particle travelling 
in a fluid. The evolution of such a particle is described by the 
Fokker-Planck equation with the diffusion tensor depending 
on the correlations of the fluctuating force in the unper- 
turbed state (e.g. Isihara 1971). These fluctuations in the 
force are generated as a consequence of the graininess of the 
background, which may be either collisional or collisionless. 

Fo r coll isionless syst ems, Chandr asekh ar & von Neu- 
mann (|1942|) , Lee (|l96c|), K andrup ( |l983| ), Bekenstein & 
Maoz ( |1 992; ) and Colpi ( ]1998; ), among others, computed the 
autocorrelation tensor in a stellar system by making the ap- 
proximation that each background particle follows a linear 
trajectory with constant velocity. The inclusio n of curve d 
orbits does not change the result appreciably (Lee 1968). 
Assuming a M axwel lian distribution of field p articl e veloci- 
ties, Kandrup (1983) and Bekenstein & Maoz (1992) derived 
a formula identical to the Chandrasekhar formula from the 
correlation tensor. 

For collisional systems, it is usually assumed that the 
mutual encounters between field particles destroy the corre- 
lation of each field particle's orbit with itself. However, from 
this it does not follow that scattering of field particles re- 
duces the correlation t ensor and so the dynamical friction. 
For instance, Ostriker (1999) considered the dissipative drag 



force experienced by a perturber of mass M p with veloc- 
ity V moving through a homogeneous gaseous medium of 
sound speed c s and unp erturbed densi ty po- Both the lin- 
earized Euler equations ( |Ostriker 1999 ) and numerical cal- 
culations (Sanchez-Salcedo & Brandenburg 1999, hereafter 
referred to as Paper I) show that the gravitational drag in 
such a fully collisional medium exceeds the value given by 
the Chandrasekhar formula, provided that the Mach num- 
ber, M = V/c B , lies in the range 1-2.5. 

In order to ease visualization it is convenient to con- 
sider the gaseous medium as a system of a large number of 

t For the drag force in the case of an absorbing perturber we 
refer to the paper by Ruffert (1996). 
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interacting field particles of identical masses m. In the hy- 
drodynamical limit the distribution of velocities is always 
Maxwellian (which is not true for collisionless systems). Let 
us assume that in such a limit the interaction between field 
particles is large enough that the effective velocity in each 
direction is close to c s for the field particles. An estimate of 
the force in gaseous media emerges immediately by extend- 
ing Chandrasekhar's results. For collisionless backgrounds, 
Chandrasekhar (1943) showed that ambient particles with 
speeds lower than the object contribute to the drag force as 
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4nG 2 M 2 m 
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dvf(v) In A + In 



V 2 
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(1) 



whereas field particles moving faster than the perturber con- 
tribute with the lower-order term 
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In the equations above f(v) dv is the number of field parti- 
cles with velocities between v and v + dv, and A = & max /&min, 
with femax the characteristic size of the medium and fe m i n = 
max{r m in,GM p /V 2 }, where r m i n is the characteristic size 
of the perturber. Now, if we blindly use the above equa- 
tions for a background of particles f(v) = NoS(v — c s ), with 
No = po/m, and keep only the dominant term for V > c s , 
we get 



F di = Ai 



4-kG 2 M; Po 
V 2 



\n(-p-)+0(ln[l-M- 2 ]), (3) 



where again M = V/c s and fe ma x ~ Vt was used. A discus- 
sion about th is dep endence of 6 ma x can be found in Ostriker 
& Davidsen (1968|). For V < c s we have 
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where Ai and A2 are dimensionless parameters of order unity. 
The first result is that the drag force is nonvanishing even for 
subsonic perturbers, as occurs in collisionless backgrounds. 
We note that the lower-order terms in Eq. (Q) are often 
ignored in the literature, which has sometimes led to wrong 
premises. For instance, Zamir (1992) ignored such lower- 
order terms and concluded that an object launched in a 
homogeneous and isotropic background with a cutoff in the 
velocity distribution from below could even be accelerated. 

We may compare the above estimates with the exact 
values obtained by c ompu ting the gravitational wake behind 
the body. Ostriker (1999) showed in linear theory that the 
drag force is given by: 



F df = 
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for M = V/c s > 1 and t > r min / (V — c s ), and 
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(6) 



for M < 1 and t > r- m i n /(c s —V). It was assumed that the 
perturber is formed at t = 0. A minimum radius r m i n was 
adopted in order to regularise the gravitational potential of 
a point mass. We see that the expressions (^) and (^) match 
Eqs (jsj) and ([]), respectively, if we take Ai = 1 and A2 = 1/2 



0. This suggests that many features of the dynamical friction 
for stellar systems, such as the role of self-gravity or the 
importance of non-local effects, may be shared by gaseous 
systems. 

In the next section we analyse the problem of orbital 
decay in a spherical and gaseous system. For stellar systems 
it was shown that Chandrasekhar's formula is a good ap- 
proximation (e.g., Lin & Tremaine 1983; Cora et al. 1997). 
It is interesting to check whether Ostriker's formula works 
also in such geometry, in analogy to stellar systems. 



2.2 Comparison with a totally absorbing 
perturber 

Considerable work has been devoted to studying the hy- 
drodynamics of the Bondi-Hoyle-Lyttleton (hereafter BHL) 
accretion problem, in which a totally absorbing body inter- 
acts gravitationally with the surrounding gaseous medium. 
In this case two different forces acting upon the accretor 
can be distinguished: the gravitational drag caused by the 
asymmetry distribution of mass, and the force associated 
with the accretion of momentum by the perturber. This 
case has been considered numerically by Ruffert (1996). Al- 
though the case of a totally absorbing accretor is different 
from the non-absorbing one considered in the present paper, 
it is interesti ng to compare the two cases. 

Ruffert (1996) showed that there is a gravitational drag 
on a totally absorbing body that saturates very quickly in 
the subsonic regime, and is roughly one order of magnitude 
smaller than in the supersonic models. This is qualitatively 
similar to the case without accretion, as described by Eqs 
(^) and (^), although the exact strength of the gravitational 
drag may be somewhat different. At first glance, the gravi- 
tational drag seems to saturate or even decline with time in 
Ruffcrt's supersonic experiments. By contrast, in both Os- 
triker's analysis and in our numerical simulations with no 
mass accretion the gravitational drag clearly increases loga- 
rithmically in time. 

Ruffert focuses mainly on the case -R so ft <S Rac, where 
_R ac = 2GM P /V 2 is the accretion radius. Whilst that case is 
relevant to BHL accretors, we are mainly interested in the 
opposite case of small accretion radii. In order to facilitate 
comparison, we will co mpar e the values of the gravitational 
drag found by Ruffert (1996) for his largest accretor models, 
EL, FL, GL and HL, where the size of the accretor is one 
accretion radius, with the values given by Eqs (^J) and ^j. 

For large accretors, the hydrodynamical force due to 
accretion acts also as a drag force (e.g. Ruffert 1996). That 
situation corresponds to geometrical accretion of mass which 
tends to suppress the flux of mass from the forward side of 
the accretor to the back side in the supersonic case. This fact 
smears out the asymmetry of the mass density distribution 
between both sides and, consequently, the gravitational drag 
is expected to be reduced. 

Ruffert (1996) gives the evolution of the gravitational 
drag force in units where J? ac = Cs = po = 1. In order to 
obtain the result in units of A-kG 2 M 2 po/ c 2 , we have to di- 
vide Ruffert's values by ivM 4 . Thus, for his model GL with 



S Note that the lower-order term of equation (Q) , which is of order 
0(ln~ 1 A), must have a coefficient 1/2 to match Eq. Q5I) . 
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M = 3 the value 80 at time t = 50i? ac /cs (see his Fig. 5b) 
corresponds to 0.31. For his model HL with A4 = 10 the 
value 800 at time t = 157?ac/c s (see his Fig. 7b) corresponds 
to 0.025. These values are somewhat smaller than those ob- 
tained by using r m i n approximately the size of the accretor in 
Eq. (^), i.e. 0.55 and 0.045, respectively. In order to trace the 
origin of these discrepancies we have checked in our records 
the drag force for M — 3 and i? a oft = Rue (a resolution of 
400 x 800 zones was used, with a mesh width at the body 
position of 1/8 accretion radius). At time t = 50-R a c/c s we 
obtained a dimensionless drag force of 0.37. The remain- 
ing 20% discrepancy may be associated with the fact that 
Ruffcrt's perturbers are absorbing. This confirms our expec- 
tations that the gravitational drag should be smaller for an 
absorbing perturber, but note that the total drag (gravita- 
tional plus momentum accretion drag) is stronger. 

Larger differences, however, are apparent for subsonic 
motions. In fact, the gravitational drag for a body travelling 
with M ach number 0.6 (model EL; see Fig. 2b in Ruffert 
(1996)) is a factor 4-5 larger than predicted by Eq. (|^). 
This may be a consequence of the fact that, in the sub- 
sonic case, accretion is more efficient in loading down mass 
from the forward side than from the back side. Therefore 
we conclude that, for large accretors, mass accretion would 
be important for dynamical friction in the subsonic regime. 
After this short digression we now return, however, to the 
case of non-absorbing perturbers with small accretion radii. 



3 SETUP OF THE MODEL 

We have performed several 3D simulations to study the dy- 
namical friction effect on a rigid perturber (or satellite) or- 
biting within a partly gaseous and initially spherical sys- 
tem, which is modelling the primary or central galaxy, with 
gas density p(r,t). The equation of state is that of a per- 
fect gas with specific heat ratio 7 = 5/3. Firstly, we will 
assume that the gas is polytropic, P oc p r , where the poly- 
tropic index V ranges between 1 and 5/3 = 7, the limits 
corresponding to the gas being isothermal and adiabatic, 
respectively. The perturber is modelled by a rigid Plum- 
mer model with the corresponding gravitational potential 
<f)p — —GMp/^Jr 2 + Rg aft , where G is the gravitational con- 
stant. In the absence of the perturber, the gas is initially 
in hydrostatic equilibrium in the overall potential given by 
4>* + cj> s , where <f> s is the potential created by the gas com- 
ponent and 0* by the rest of the mass. The potential tj>* is 
assumed to be fixed in time, and the centre of the primary 
fixed in space. Self-gravity of the gas is ignored which seems 
to be a reasonable approximation to reality (e.g. Prugniel & 
Combes 1992). For simplicity we assume that the combined 
potential 0g = <t>* + 4>s is also given by a Plummer potential 
with softening radius _Rg and total mass Mg- As a conse- 
quence, a polytropic system of pure gas must have polytropic 
index P — 1.2 in order to have a self-gravitating model (e.g. 
Binney & Tremaine 1987). For F > 1.2 the models above are 
not self-consistent at large radii because the gravitational 
force due to the gas increases faster than |V^g|- However, 
we use simple initial condition since our aim is to test the 
applicability of various fit formulae which should be able 
to predict the evolution of the satellite regardless of the full 
consistency of the background system. For completeness, the 



case of 0g being a King model instead of a Plummer model 
is also explored in order to mimic the inner portions of a 



self-gravitating sphere; see Section 4.5 



Using an explicit cartesian code that is sixth order in 
space and third order in time, we solve numerically the con- 
tinuity and Euler equations which, for a polytropic gas, are 

dp 
dt 

dv 

dt 

where h is related to the specific enthalpy and is defined by 



+ V ■ (pv) -- 
+ (v ■ V)w 



0. 



-V(/l + 4>g + 4> P ) + viscous force, 
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r-i 



1 



Csoln(p/po) 



if T / 1 
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(9) 



where c s o = y^Po/po, Po and po are the reference values of 
pressure and density, respectively, and 



GMt, 



^(r-R p (t))* + Rl it 



(10) 



is the gravitational potential generated by the perturber at 
the position R p (t). Like in Paper I, an artificial viscosity 
has been introduced in the momentum equations. Further 
details about the numerical scheme adopted are given in the 
appendix. 

We will also investigate the more general case in which 
Eq. (bh is replaced by 



V(0 G + < 



viscous force, (11) 



(12) 



dv . VP 

— + (v-V)v = 

at p 

and the specific entropy, s, evolves according to 
ds 

— + (v ■ V)s = viscous heating. 

The perturber, which is introduced instantaneously at 
t = 0, moves in the (x, j/)-plane and evolves according to the 
equation 

Mp ^W~ = ~ m p v<?!>g + Fdf ' ( 13 ) 

with the dynamical friction force, Fdf, given by 



F df = / 5pV<t> p (fr, 



(14) 



where Sp = p(r, i) — p(r, 0). In all the simulations presented 
here V ■ Fdi < and Fdi ■ e T < at any time of the run, 
where e r is the unit vector in the radial direction. 

A symmetry condition is applied at z = 0. Apart from 
that the domain is a rectangular box with open boundary 
conditions. The size of the box must be taken large enough to 
ensure that the density perturbations that have propagated 
outside the domain do not contribute significantly to the 
friction integral (|l4|). To improve matters, additional spa- 
tial extend has been gained by implementing a non-uniform 
mesh (see appendix for details) . Most of the simulations were 
carried out with the grid represented in Fig. ^. 

We choose units such that G = Rg = p(0, 0) = 1, where 
r = corresponds to the centre of the primary and t — 
is the beginning of the simulation. The parameters to be 
specified in the polytropic case are P, the isothermal sound 
speed Cso, the mass of the central galaxy Mg, the mass M p 
and softening radius of the satellite, and its initial distance 
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Figure 1. Coordinates for meshpoints in the non- uniform grid 
with a transformation as given by Eq. ( |Al| ), The dotted line cor- 
responds to the z-coordinate and the solid line for the x and 
y-coordinates. 



to the centre, and velocity. We do not try to simulate the 
complex case in which the accretion radius is larger than 
the softening radius. That would require a detailed stability 
analysis of the flow. Objects with softening radius larger 
than their accretion radius, like galaxies in clusters, might 
only suffer from fli p-flop instabilities if a large amount of 
gas is replenished (Balsara, Livio & O'Dea 1994). It is still 



unclear whether this instability occurs in 3D or not. 

Two requirements are imposed in the initial distance of 
the satellite. Firstly, the initial distance is required to be 
a few galactic core radii (-Rg) and, secondly, we will con- 
sider orbits such that i> c (-Rapo)/c s (i? a po) > 1, where 7? ap o is 
the apocentric distance of the perturber's orbit in the ab- 
sence of drag, v c (r) = */ rd<f>G / dr the circular velocity and 
Cg = TP/ p for a polytrope. The first condition is necessary 
in order to have the complete history of the evolution of 
the merging of an 'external' satellite. The second condition 
comes from the virial theorem as follows. Applying the virial 
theorem for a spherical system with no fluid motions we get 



V + 3(r-l)W- 



= 0, 



(15) 



where V = — J"' s pv\ d 3 r and U is the total internal energy 
within the sphere of radius rs. We immediately obtain that 



(16) 



(17) 



He r dp 
cf p dr 

for a polytrope, and 

Vc r dp 

cf 7p dr 

for an isothermal background, assuming that perturbations 
are adiabatic. Since the density should decay like 1/r 2 or 
faster (r~ n with n > 2) in the outer radii, the Mach number 
satisfies M = v c /c s > in the outer parts. Therefore, 

for a monatomic gas, circular orbits in the outer parts are 
always supersonic. 



Figure 4. Different components and strength of F^f versus time 
for the reference model, F^f (solid line), azimuthal component 
of F d f (dotted line) and radial component of f d f (dot-dashed 
line) for a freely decaying perturber. We have taken the absolute 
values but it is worthwhile to recall that both F^f ■ V and F^t ■ e r 
have negative values. The temporal evolution of _F<jf r (triple dot- 
dashed line) and Fdf ib (dashed line) for a perturber forced to 
rotate in a circular orbit with R p (0) = 4 are also plotted. 



4 RESULTS 

4.1 The reference model 

In this subsection the evolution of the perturber is discussed 
in a reference model. Variations from this model are consid- 
ered in subsequent subsections. Our reference model uses a 
resolution of 95 x 95 x 40 meshpoints, and a non-uniform grid 
as shown in Fig. |l|. The gas is assumed to be a polytrope 
with r = 1. The satellite has a mass M p = 0.3, softening 
radius -R so ft = 0.5, and an initially circular orbit of radius 
R P (0) = 4. We have chosen Mq = 10 and c s = 1.3, so the 
mass of gas within R P (0) is approximately Mg Q Thus, the 
satellite mass is ~ 1/30 of the galaxy mass, and the soften- 
ing radius is three times the Bondi radius. Notice that the 
perturber will move supersonically until a radius of approx- 
imately 0.48. 

The satellite's angular momentum per unit mass, J p , 
and orbital radius are shown in Fig. |as a function of time. 
We see that the angular momentum decreases almost lin- 
early with time (solid line). A snapshot of the density en- 
hancement and the velocity field at the plane z — and 
t = 26 is shown in Fig. |^. The strength of the radial and 
azimuthal components of Fat, in units of 4irG 2 M^po/ 'c% 
(po = p(r, 0) hereafter) evaluated at the instantaneous po- 
sition of the satellite (the 'dimensionless force' hereafter), 
are presented in Fig. [| for the standard model. As expected, 
the dimensionless force saturates since the direction of V 
changes and, consequently, the wake behind the body effec- 
tively 'restarts'. The force shows some oscillations caused by 
the combination of the interaction of the satellite with its 



" Typical parameters for galaxies are Rq 
g/cm 3 , so the mass unit is 10 10 M©. 



10 kpc, po = irr 
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Figure 2. Evolution of the satellite in the reference model. The time unit is (Gp(O,0)) In panel (a) is plotted the orbit in the 

(x, i/)-plane. The angular momentum (per unit mass) and orbital radius of the perturber are shown as solid lines in panels (b) and (c), 
respectively, together with the predictions of the Local Approximation Prescription (LAP; dotted lines). Symbols mark the time when 
the orbital phase is zero, so the time between two consecutive symbols is one orbital period. 



own wake, together with small epicyclic motions of the satel- 
lite. For a perturber forced to rotate in its originally circular 
orbit at constant velocity, the azimuthal force is constant 
with time (dashed line in Fig. ^|), reaching a plateau after 
half an orbital period. Almost no oscillations in the force are 
seen in this case. 

Let us compare the decay rate with the prediction of 
formulas (||) and (|^). Denoting the azimuthal component of 
the force by Fdi,ip, this expression yields for R p > 0.48 (i.e. 
in the supersonic regime): 



dJ p 
~~dt 



Fdf^R P _ 4ttG 2 M p poR p 



Vc 



(18) 



where v c is the circular velocity, and 7 a< ji = In A + 
| In (l — .M~ 2 ) a dimensionless quantity. From J p = v c R P 
we have 



dRp 
~~dT 



2-1 



R r , 



2RI + 1 



dJ p 
dt 



(19) 



It is noted that now po and M depend on the position. In 
order to integrate Eq. ([ll]), the only remaining unknown is 
In A. Since In A is no longer a linear function of time, we 
assume that 7 a di takes the following values for Ai > 1 



where 



ln(-^) + |ln5 A1 



S M = l- M~, t c 



if r <t <t c 
if t <t <t c 
if r <t c <t 



(20) 



(21) 



V c — C s V c + C s 

with 13 being a free parameter which is chosen to match the 



values of J p (t) and R P (t) obtained numerically. The mini- 
mum radius is taken to be r m i n = 2.257? ao ft, as was inferred 
for homogeneous media (Paper I) provided that the soft- 
ening radius is larger or equal to the accretion radius. For 
t < t < t c a linear interpolation has been adopted. For times 
such that t c < t, it is difficult to suggest any prescription 
from the linear theory. We therefore assume that the value 
remains approximately unchanged compared to the value it 
has got just before coming into that interval. However, this 
assumption is not completely satisfactory for bodies moving 
initially at Mach numbers close to 1. Therefore we addition- 
ally assume that 



iadi = — 
T 



hi 



( — ) 

\M - 1J 



+ - In Sm 



(22) 



if t < t and Ai > 1, regardless of the condition r < t c . 

For M < 1 and a homogeneous background the lin- 
ear theory is not able to capture the temporal evolution of 
the drag force but just the asymptotic values; see Eq. (^). 
As a consequence, for a body moving on a curved orbit, no 
reliable estimate of the drag force can be proposed straight- 
forward. In the following we shall refer to Eq. ^ together 
with Eqs (^),(^2|) as the Local Approximation Prescription 
(LAP). It is 'local' in the sense that the force depends only 
on quantities at the instantaneous position of the satellite 
even though, of course, the wake extends far away beyond 
the satellite. At this stage it becomes clear that the LAP 
may fail at least for M close to or less than unity. 

The results of the integration for the reference model 
are plotted as dotted lines in Fig. ^|for j3 = 3.1. We see that 
it agrees closely with the numerical values. Part of the suc- 
cess of Eqs (^)-(^) resides in the fact that the condition 
t < t c is satisfied very soon, at R p ~ 3.5, for the present 
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Figure 3. Grey-scale plot of the density enhancement In (p(r, t)/p(r, 0)), together with velocity vectors, at z = and t = 26, for the 
reference model with M p = 0.3 and i? so f t = 0.5. 



case. However, it may be somewhat different for perturbers 
moving with Mach numbers close to 1 for a signifi cant in- 
terval of time. This case will be considered in Sect. 



4.3 



We ran a model which was identical to the reference 
model except that the perturber mass was increased to 



0.5. The results of the simulation are plotted in 



Fig. ^a, together with LAP curves of different (3. There is 
very good agreement with the LAP for /3 = 3.1. In general, 



it is found that a change of a factor three in the mass of the 
perturber produces a slight variation of about five per cent 
in the azimuthal component of the dimensionless drag force. 
Fig. ^a also shows J p (t) for the same perturber but initially 
at orbital radius R p = 6.0. Once the angular momentum has 
decayed to 6.0 the further evolution is quite similar to the 
evolution of J p starting at 6.0. 



A new run was done for Af D 



0.5 and R s , 



1. The 
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/3cir = 



"min 1 + M 1 

r p {i-M- 2 y/ 2 



exp (/adi) 



(23) 



Figure 5. (a) Temporal evolution of the angular momentum (per 
unit mass) for the isothermal model with M p = 0.5, i? so ft = 0.5 
and -Rp(O) = 4 (solid line) and for an identical perturber but 
initially at R p (0) = 6 (dashed line). Panel (b) is like panel (a) 
but now iJ so f t = 1 and R p (0) = 4. The predictions by LAP 
with different /3 are also plotted for the perturbers initially at 
R p (0) = 4. LAP overestimates the drag when the distance of 
perturber to the centre becomes comparable to R so f t - 



evolution of J p is plotted in Fig gb. Again the LAP using 
(5 = 3.1 reproduces reasonably well its evolution. 

For simulations with resolutions 120 x 120 x 54 the 
change in the drag force is less than 0.5 per cent. 

A measure of the robustness of the fitting formulae is 
given by the dispersion of (3 for forced circular orbits of 
different radii, which will be referred to as /3 i r . Using for 
/3 c ir the expression 



with Tadi being the asymptotic value of the azimuthal com- 
ponent of the force obtained numerically, we have /3 c ir = 
2.6, 2.9, 3.2 and 3.5 for circular orbits of radii R p = 4, 3, 2 
and 1.5, respectively. The rotation velocities were taken ac- 
cording to the rotation curve of the Plummer model. If the 
perturber moves on a circular orbit of radius R p = 4 but 
M — 1.61, ,0cir is found to be 3.6 (see Table 1). These re- 
sults suggest that /3 c i r depends on M, but could also contain 
information about the global density profile. In the inner re- 
gions /3 c ir should also contain the effect of the excitation 
of gravity waves. Strong resonances occur when the local 
Brunt- Vaisala frequency matches the orbital frequency of 
the perturber. The existence of resonances depends strongly 
on the temperature profile in the inner regions and require 
in general a large-scale entropy gradient in that region (Bal- 
bus & Soker 1990; Lufkin et al. 1993). For satellite galaxies 
the amplification of gravity waves, if any, may only take 
place in the very inner part of the parent galaxy. However, 
other effects, such as mass stripping and tidal deformation 
of the satellite galaxy as well as friction with the luminous 
and collisionless parts of the galaxy, may play an important 
role as far as dynamical evolution of the satellite galaxy is 
concerned. 

For fixed rotation curve, the parameter /? may have 
some dependence on other variables of the system such as 
r. This question will be considered in the next few sections. 



4.2 Testing the accuracy. Constant-velocity 

perturbers in a homogeneous background 

A measure of the accuracy of our 3D simulations is given 
by comparing the force of dynamical friction experienced 
by a perturber moving on a straight-line orbit with con- 
stant velocity V in a uniform and isothermal medium, with 
the axisymmetric simulations described in Paper I, which 
had higher resolution. In Fig. |^ we plot the drag force ex- 
perienced by a body that moves in the i-direction using 
the same resolution as that in our reference model, together 
with the values obtained in high-resolution simulations. The 
initial position of the body is (—4, —4, 0) and it travels at 
constant velocity with Mach number 0.75 (bottom lines) and 
1.12 (upper lines). In the supersonic experiments the ratio 
between Rsoft and the accretion radius defined as 2GM P /V 2 , 
was 1.5, with a softening radius of 0.5. As long as the body is 
within the region |a;| < 6 the accuracy of the force is better 
than five per cent. 



4.3 Poly tropic gas 

Let us now consider the same perturber with i? S oft =0.5 and 
mass M p that is forced to move on a circular orbit in a poly- 
tropic gas and compute the value of /3 c i r that matches the re- 
sults. We do this for models which have different poly tropic 
indices F. The rotation curve is assumed to be the same 
as in the reference model. Rotation curve and sound speed 
profiles for models with different pairs (r, c a o) are drawn in 
Fig. [?]. Since they have different polytropic indices we will 
specify each of these models by just giving their polytropic 
index. For the model with polytropic index 5/3, the Mach 
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Figure 6. Dimensionless drag force experienced by a body mov- 
ing according to R p = (—4 + Vt, —4, 0) in a unperturbed back- 
ground of constant density. The solid line shows the numerical 
results for high resolution axisymmetric simulations whereas the 
dotted line is for 3D simulations with resolution as in the reference 
model. 







I \ 7 p =5/3 

























Figure 7. Rotation curve for the Plummer model (solid line), and 
the sound speed profiles for T = 5/3 (dot-dashed line), T = 1.4 
(triple dot-dashed line), and T = 1.2 (dashed line). 



as defined in Eq. ( pj| ) for models which have dif- 
ferent r. The body has a mass of 0.3, softening radius 0.5, and 
is forced to move in a circular orbit of radius R p with a Mach 
number which may not correspond to the rotation curve of the 
Plummer model. If /3 c ; r are displayed in the same row for both 
Rp = 4 and R p = 2, it means that they correspond to the rotation 
curve given by the model. 



r 


M 


/3 cir 


M 


/3 cir 




Rp =4 


Rp =4 


R P = 2 


Rp = 2 


1.0 


1.16 


2.6 


1.45 


3.2 


1.0 


1.61 


3.6 






1.2 


1.61 


2.8 


1.67 


3.1 


1.2 


1.27 


2.1 






1.2 






1.41 


2.7 


1.4 


1.27 


2.2 


1.41 


2.7 


1.4 


1.61 


3.0 






1.4 






1.67 


3.0 


5/3 


1.12 


2.4 


1.09 





bit. The case of F = 5/3 is particularly interesting because 
then t c < t for a remarkably broad range of R p . The j3- 
parameter seems to be rather arbitrary in that range. In 
Fig. |^ the evolution of J p and R p is plotted for M p = 0.3 
(upper panel) and M p = 0.1 (bottom panel). In both sim- 
ulations r = 5/3. The prediction of the LAP is plotted as 
dotted lines for j3 = 2.1 in both cases. In the first case the 
decay is so fast that the LAP might no longer be a good ap- 
proximation (upper panel). However, even for slower decays 
(the case of Mp = 0.1) the LAP is not able to give reliable 
values of the force in the subsonic regime. In fact, Eq. Q 
overestimates significantly the friction, mainly because of 
two reasons; firstly, the wake must restart as a consequence 
of the curvature of the orbit and, secondly, the orbital radius 
and size of the object become comparable. 

Since for the model with r = 5/3 the typical velocities 
v c /cs lie in the range 1.0-1.2 at the radius interval € [2,4], 
one would expect, according to Eq. (^), a larger value of 
the dimensionless drag force in this case compared to the 
isothermal gas model in which the body moves with higher 
values of M . A comparison between those forces for T = 5/3 
and r = 1 is given in Fig. []. Unexpectedly, the dimensionless 
drag force is stronger in the run corresponding to T — 1. 

Finally, the evolution of J p of a perturber of mass 0.3 
is drawn in Fig. ^ for V — 1.4 as well as for T = 1.2. The 
/3-parameter was 2.2 and 3.1, respectively. Clearly, the LAP 
becomes inaccurate as soon as the body becomes transonic. 



number of the perturber at R p = 4 is approximately the 
same as in the reference model. The corresponding values of 
/3 c ir for each value of r are given in Table 1, except for the 
case F = 5/3 at R p = 2 because for that orbit t c <C T. Table 
1 suggests a remarkable dependence of /3 c i r on M, as well 
as on P. This reflects the fact that the asymptotic value of 
the azimuthal component of the dimensionless force is not 
sensitive to the Mach number whenever M = 1.1—1.6. Given 
P and R p , the variation of the force in that interval of Ai is 
less than five per cent. 

For sinking orbits the value of may differ somewhat 
from the median of /3 c ir because the orbit of the body loses 
correlation with its wake faster than in a fixed circular or- 



4.4 Non-circular orbits in the Plummer model 

The dependence of the orbital sinking times on the eccen- 
tricity of the orbit is of interest for the statistical analysis of 
merging rates of substructure in the cosmological scenario 
(e.g. Lacey & Cole 1993; Colpi et al. 1999), and for the the- 
oretical eccentricit y distributions of globular clusters and 
galactic satellites ( van den Bosch et al. 1999[ ). On the one 



hand, the hierarchical clustering model predicts that most of 
the satellite's orbits have eccentricities between 0.6 and 0.8 
(Ghigna et al. 1998). On the other hand, the median value 
of the eccentricity of an isotropic distribution is typically 0.6 
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Figure 8. In the left panels it is shown the evolution of the 
angular momentum (per unit mass) for the model with T = 5/3 
and M p = 0.3 (solid line of upper panel) and M p = 0.1 (bottom 
panel). The corresponding evolution of the radial distance of the 
perturber is plotted as solid lines in the right panels. Dotted lines 
show the predictions of the LAP using j3 = 2.1. LAP fails as soon 
as the motion becomes slightly supersonic. 
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Figure 9. Solid lines show the value of the radial component of 
the force F^i for the T = 1 model with M p = 0.3 (at the top) 
and for T = 5/3 and M p = 0.1 (at the bottom). The azimuthal 
components are plotted with dashed-lines, whilst the radial ones 
in dotted lines. 




Figure 10. Decay of J p of a perturber of mass M p = 0.3 initially 
in circular orbit, for different values of T. The dotted lines indicate 
the prediction of the LAP. 
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Figure 11. Panel (a) shows the orbital angular momentum in 
units of the initial value for orbits having initially the same energy 
but different eccentricities in the Plummer model with T = 1. 
M p = 0.4 and R so f t = 0.5 were taken. In panels (b), (c) and 
(d) the corresponding orbits in the plane (x, y) are presented for 
initial circularities 0.75, 0.6 and 0.5, respectively. 



( van den Bosch et al. 1999 ) . Here we explore the dependence 
on eccentricity for dynamical friction in a gaseous sphere in 
the Plummer potential. The non-singular isothermal sphere 
is considered in the next subsection. 

We shall discuss the sinking rate in terms of the an- 
gular momentum half-time, ii/2, the time after which the 
perturber has lost half of its initial orbital angular momen- 
tum. Fig. [ll] shows the dynamical evolution of the satellite 
on bound orbits having initially equal energies but different 
eccentricities (J p is in units of the initial value). Here we 
have used the polytropic model with F = 1. A good fit to 
the numerical results is oc e~ 0,4 , where e is the initial 
circularity, e = J P (E) / J p]c i T (E), the ratio between the or- 
bital angular momentum and that of the circular orbit hav- 
ing the same energy E. The dynamical friction time scale is 
found to increase with increasing eccentricity. This finding 
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makes a clear distinction between dynamical friction in stel- 
lar and gaseous systems. In fact, the dynamical friction time 
scale has been obtained by direct N-body simulations and in 
the linear response theory suggesting t\ii oc e ' 4 for bodies 
embedded in a truncated non-singular isothermal sphere of 
collisionless matter (van den Bosch et al. 1999; Colpi et al. 
1999). A notorious difference between stellar and gaseous 
backgrounds is that the gaseous force is strongly suppressed 
when the perturber falls in the subsonic regime along its 
non-circular orbit. In order to discern how much the above 
result depends on the model, we will consider in the next 
subsection the angular momentum half-time for bodies or- 
biting in the non-singular isothermal sphere. 

4.5 Sinking perturbers in the non-singular 
isothermal sphere 

The flat behaviour of the rotation curves of spiral galax- 
ies strongly suggests the existence of isothermal dark halos 
around them. That is the reason why most of the studies 
about dynamical friction assume an isothermal sphere as 
the standard background. In the case of a self-gravitating 
sphere of gas, the non-singular isothermal sphere is also a 
natural choice. For instance, we may consider the dynamical 
friction of condensed objects in early stages of a star cluster 
or a galaxy, in which most of the mass of the system is gas. 
The gas component is expected to follow the non-singular 
isothermal sphere out to great radii. 

We use the King approximation to the inner portions of 
an isothermal sphere, with softening radius Rg and central 
density pg(0), i.e. 



1 .o r^> 



-47rGp G (0)i?5-ln 

x 



(24) 



with x = r/Ra (Binney & Tremaine 1987). Our units are 
still G = Rg — p(0, 0) = 1. Even though the rotation curve 
adopted is very similar to that of the Plummer model, the 
detailed evolution, such as the dependence of the decay rate 
on eccentricity or the level of circularization of the orbits, 
may be somewhat different. For the sake of brevity, we only 
present the results for the case where the Euler and conti- 
nuity equations were solved together with the entropy equa- 
tion. 

Fig. [l^ shows the decay of a perturber with M p — 0.3 
and -Rsoft = 0.5 in a King model with pa — 1, i.e. pure 
gaseous sphere, for different circularities. There are no ap- 
preciable variations of the decay time for different initial cir- 
cularities. In order to gain deeper insight into the differences 
between the dynamical friction timescale in an isothermal 
sphere of gas instead of collisionless matter, we have com- 
puted the circularity at the times when the orbital phase is 
zero (Fig. |l3"| ) . The circularity does not change at all for ini- 
tially circular orbits. For eccentric orbits, deviations in the 
circularity occur but there is no net generation. The lack of 
dependence of the decay times on circularity, together with 
the absence of any significant amount of circularization, lead 
us to conclude that dynamical friction in a gaseous isother- 
mal sphere is not able to produce changes in the distribution 
of orbital eccentricities. Mass stripping which could lead to 
circularization has been ignored. 
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Figure 12. The orbital angular momentum in units of the initial 
value for orbits having initially the same energy but different 
circularities in the King model. 
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Figure 13. Cir cula rity as a function of time for the same situ- 
ation as in Fig. rL2l The circularities are computed at the time 
when the orbital phase is zero. 



5 DISCUSSION AND CONCLUSIONS 

In this paper we have presented simulations of the decay of 
a rigid perturber in a gaseous sphere. As a first approxima- 
tion, both the mass stripping and the barycentric motion of 
the primary were neglected. In addition, the softening radii 
of the perturbers were taken a few times larger than the 
accretion radius. Simulations of the evolution of very com- 
pact bodies ( JZ so f t <C -Rac), such as black holes in a spherical 
background are numerically very expensive. 

The Linear Approximation Prescription is able to ex- 
plain successfully the evolution of the perturber provided 
the (i) decay is not too fast, (ii) the motion is supersonic and 
(iii) the body is not too close to the center to avoid that R p 
becomes comparable to R so ft- In our simulations, the angu- 
lar momentum of the perturber decays almost linearly with 
time. Generally speaking, the associated Coulomb logarithm 
should depend on the distance of the perturber to the centre 
and on its Mach number. However, the evolution of a free 
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perturber is relatively well described by just a linear depen- 
dence of the maximum impact parameter of the Coulomb 
logarithm on radius, except very close to the centre where 
the time averaged Mach number is less than unity. 

In Chandrasekhar's formalism the Coulomb logarithm 
is often approximated by the logarithm of the ratio of the 
maximum to minimum impact parameters of the perturber. 
In a singular isothermal sphere of collisionless matter, the 
time for an approximately circular orbit to reach the center 
is 



1.2- 



GM P In §■ 



(25) 



where fo ma x is roughly the virial radius of the primary and 
fomin is of the order of the softening radius of the perturber 
(Binney & Tremaine 1987). We have confirmed numerically 



Ostrikcr's suggestion that the Coulomb logarithm should be 
replaced by 



In- 



1 



0.428 



■In 



R P (t) 

Rsoft 



(26) 



for the gaseous case, with r\ ~ 0.35, and independent of the 
eccentricity. For very light perturbers in near-circular orbits 
(objects with mass fractions smaller than one per cent), the 
time scale for dynamical friction is a factor 2 smaller in the 
case of a gaseous sphere than in corresponding stellar sys- 
tems (e.g. globular clusters in the Galactic halo). This factor 
is somewhat smaller when _R P becomes comparable to a few 
.Rsoft • The weak dependence of the decay times on eccentric- 
ity suggests that more eccentric orbits of merging satellites 
should not decay more rapidly in the halo during the epoch 
in which matter is mainly gas. Thus, they do not touch the 
disc at an earlier time than less eccentric orbits. Our simu- 
lations give support to the idea that dynamical friction does 
not produce changes in the distribution of orbital eccentric- 
ities. 

The analysis must be modified if one wishes to account 
for the existence of clouds in a non-uniform medium. So 
far, in a variety of systems it is assumed that the per- 
turber has condensed from fragmentation of the gas; thus, 
the gaseous component must be thought of as being dis- 
tributed in clumps or clouds even more massive than the 
'perturbers'. In certain situations, compact perturbers may 
not be dragged but heated (e.g. Gorti & Bhall 1996). For 
this reason one may be pessimistic towards a recent sug- 
gestion by Ostriker (1999) that, because of the contribution 
of the gaseous dynamical friction force, young stellar clus- 
ters should appear significantly more relaxed than usually 
expected. 
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APPENDIX A: DESCRIPTION OF THE CODE 

We solve Eqs (f^)~(|l^) on a nonuniform cartesian mesh, 
(x,y,z). This is accomplished using a coordinate transfor- 
mation to a uniform mesh, (£,r/, £), with 

*=^ ) = wfor (A1) 

where is a length parameter, and similarly for y = y(r]) 
and z — z(Q. The transformed equations are solved using 
sixth order centred finite differences, 
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Table Al. Coefficients for the derivative formulae 
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c-p x 60 
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cf ] x 180 


-490 


+270 


-27 


+2 



[wjr^Pf^ (A2) 

for th e nt h ^-derivative. The coefficients c'. are given in 



Al 



Table 
analogous to Eq. A.2 



The expressions for the n and £ derivatives are 



The corresponding x-derivative of a function f(£(x)) is 
obtained using the chain rule, so 

df _ d£df ± 

dx dxdi x" 1 ' 

where primes denote ^-derivatives. For the second derivative 
we have 

8 2 f _ u"x' - u'x" 

dx' 2 x' 3 ( ' 

Again, the expressions for the y and z derivatives are analo- 
gous. Since the scheme is accurate to A£ 6 , i-derivatives are 
accurate to Aa; 6 (£), which varies of course across the mesh. 

The third order Runge-Kutta scheme can be written in 
three steps (Rogallo 1981): 

1st step: / = / + JiAtf, g = f + CiAtf, 

2nd step: / = g + j 2 Atf, g = f + £ 2 A*/, (A5) 

3rd step: / = g + 73 At/, 

where 

8 5 3 , 17 . 5 

71= 15' 72 =12' 73 = 4' Cl = "60' ^--U^ 

Here, / and g always refer to the current value (so the same 
space in memory can be used), but / is evaluated only once 
at the beginning of each of the three steps. The length of the 
time step must always be a certain fraction of the Courant- 
Friedrich-Levy condition, i.e. At = fccFL Aa;/C/ ma x, where 
fccFL < 1 and f/max is the maximum transport speed in the 
system. 

The nonuniform mesh allows us to move the boundaries 
far away, so the precise location of the boundaries should 
not matter. In all cases we have used open boundary con- 
ditions by calculating derivatives on the boundaries using a 
one-sided difference formula accurate to second order. This 
condition proved very robust and satisfied our demands. 

Equations (p])-(|l2|) are solved in non-conservative form. 
This is sufficiently accurate because of the use of high order 
finite differences and because the solutions presented in this 
paper are sufficiently well resolved. 
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